classdef GridDatasetClass < handle
    %GRIDDATASETCLASS 此处显示有关此类的摘要
    %   此处显示详细说明
    
    properties
        filePath
        Allfilelist
    end
    
    methods
        function obj = GridDatasetClass(filePath)
            obj.filePath = filePath;
            yearlist = dir(obj.filePath);
            yearlist = yearlist(3:end);
            Allfilelist = [];
            for i = 1:length(yearlist)
                filelistbyyear = dir([filePath '\' yearlist(i).name]);
                filelistbyyear = filelistbyyear(3:end);
                Allfilelist = [Allfilelist;filelistbyyear];
            end
            obj.Allfilelist = Allfilelist;
        end




        % 对原生nc文件,求rho,dh并联合pt,s保存到mat
        function [lon, lat, depth, pt, s, rho, dh] = rho_dh_measure(obj,d)
            filelist = obj.Allfilelist;
            file = [filelist(1).folder '\' filelist(1).name];
            out = GridDatasetClass.ncFileRead(file, 'longitude','latitude','depth');
            lon = out.longitude;
            lat = out.latitude;
            depth = out.depth;
            for filei = d
            % parfor filei = 523:length(filelist)
                % filei
                file = [filelist(filei).folder '\' filelist(filei).name];
                if ~exist(['./' filelist(filei).folder(end-3:end) '/'], 'dir')
                    mkdir(['./' filelist(filei).folder(end-3:end) '/']);
                end
                out = GridDatasetClass.ncFileRead(file,'thetao','so');
                [pt,s,rho,dh] = GridDatasetClass.potiential_temp_measure_from_temp_sail(lon, lat, depth, out.thetao, out.so);
                savefilename = ['E:/code/GOPRncFile2/' filelist(filei).folder(end-3:end) '/' filelist(filei).name(1:end-3) '.mat'];
                GridDatasetClass.parsave(savefilename, lon, lat, depth, pt, s, rho, dh);
            end
        end

        function ssta_measure_by_each_year_month(obj)
               
            %创建分组每组个数索引 
            MonthGroupCount = zeros(12);

            filelist = obj.Allfilelist;
            for i = 1:length(filelist)
                % 提取文件名字 格式为 yyyymmdd.nc 如20080101.nc
                filename = filelist(i).name;
                format = 'yyyyMMdd';
                date = datetime(filename(1:end-4),'InputFormat',format);
                Month = date.Month;
                % 对该i索引文件对应的Group数量加1，并添加i进Idx矩阵
                MonthGroupCount(Month) = MonthGroupCount(Month) + 1;
                MonthGroupIdxMatrix(Month, MonthGroupCount(Month)) = i;
            end
            [m,~] = size(MonthGroupIdxMatrix);
            for i = 1:m
                groupfilelistidx = MonthGroupIdxMatrix(i,:);
                groupfilelistidx(groupfilelistidx==0) = [];
                groupfilelist = filelist(groupfilelistidx);
                GridDatasetClass.group_measure_and_save(groupfilelist);
            end

        end

        function ssta_measure_by_each_year_day(obj) 
            %创建分组每组个数索引 
            MonthGroupCount = zeros(366);

            filelist = obj.Allfilelist;
            for i = 1:length(filelist)
                % 提取文件名字 格式为 yyyymmdd.nc 如20080101.nc
                filename = filelist(i).name;
                format = 'yyyyMMdd';
                date = datetime(filename(1:end-4),'InputFormat',format);
                Month = date.Month;
                Day = date.Day;
                % 随便取一个润年
                d = days(datetime(2000,Month,Day) - datetime(2000,1,1)) + 1;
                % 对该i索引文件对应的Group数量加1，并添加i进Idx矩阵
                MonthGroupCount(d) = MonthGroupCount(d) + 1;
                MonthGroupIdxMatrix(d, MonthGroupCount(d)) = i;
            end
            [m,~] = size(MonthGroupIdxMatrix);
            for i = 1:m
                groupfilelistidx = MonthGroupIdxMatrix(i,:);
                groupfilelistidx(groupfilelistidx==0) = [];
                groupfilelist = filelist(groupfilelistidx);
                GridDatasetClass.group_measure_and_save(groupfilelist);
            end

        end


        function ssta_measure_by_season(obj, varargin)
            filelist = obj.Allfilelist;
            % filelist = filelist(9983:end);
            % 设置分组起始锚点
            begintensor = 1;

            % 设置滑动锚点
            datetensor = begintensor;

            % % 设置分组结束锚点
            % endtensor = 0;

            % 添加分组矩阵
            idx_groupArray = find(strcmp(varargin,'groupArray'));
            if ~isempty(idx_groupArray)
                seasonArray = varargin{idx_groupArray + 1};
            else
                seasonArray = [12 1 2;3 4 5;6 7 8;9 10 11];
            end
            
            [m,~] = size(seasonArray);
            % 添加分组索引
            seasonid = 1;

            % 通过索引取出当前分组季节
            season = seasonArray(seasonid,:);
            
            % 通过锚点批量处理filelist
            while 1
                % 提取文件名字 格式为 yyyymmdd.nc 如20080101.nc
                filename = filelist(datetensor).name;
                format = 'yyyyMMdd';
                date = datetime(filename(1:end-4),'InputFormat',format);


                % 如果月份属于该季节继续，否则分组完成
                if ~ismember(date.Month, season) 
                    % 当分组完成时,记录结束锚点
                    endtensor = datetensor - 1;
                    % 判断endtensor是否大于等于begintensor，若小于则无分组
                    if endtensor >= begintensor
                        %取出分组
                        groupfilelist = filelist(begintensor:endtensor);
                        GridDatasetClass.group_measure_and_save(groupfilelist)

                        %重置锚点
                        begintensor = endtensor + 1;
                        datetensor = begintensor;

                        %季节递进
                        seasonid = seasonid + 1;
                        if seasonid > m
                            seasonid = seasonid - m;
                        end
                        season = seasonArray(seasonid,:);
                    else
                        seasonid = seasonid + 1;
                        if seasonid > m
                            seasonid = seasonid - m;
                        end
                        season = seasonArray(seasonid,:);

                    end
                else
                    datetensor = datetensor + 1;
                end

                % 判断越界强制结束
                if datetensor > length(filelist)
                    endtensor = datetensor - days(1);
                    groupfilelist = filelist(begintensor, endtensor);
                    GridDatasetClass.group_measure_and_save(groupfilelist);
                    break;
                end
            end


        end

        % function mean_pt_measure_by_day(obj)
        %     filelist = dir(obj.filePath);
        %     filelist = filelist(3:end);
        % 
        %     Year = 2000;%找一个闰年
        %     startDate = datetime(Year,1,1);
        %     tensorDate = startDate;
        %     endDate = datetime(Year,12,31);
        %     GDC = GridDatasetClass(obj.filePath);
        % 
        %     while tensorDate <= endDate
        %         Sum_s = 0;
        %         Sum_s_effictive = 0;
        %         Sum_t = 0;
        %         Sum_t_effictive = 0;
        % 
        %         for i = 1:length(filelist)
        %             i
        %             if ~exist(['./' filelist(i).name '/'], 'dir')
        %                 mkdir(['./' filelist(i).name '/']);
        %             end
        %             date = datetime(str2double(filelist(i).name),tensorDate.Month,tensorDate.Day);
        %             try
        %                 GDC.GridDataSelect(date);
        %             catch
        %                 warning([datestr(date,'yyyymmdd') '未找到文件']);
        %                 continue;
        %             end
        %             out = GDC.ncFileRead('longitude','latitude','thetao','depth','so');
        % 
        % 
        %             %%%%%%%%%%%%%%%%%% 根据实际文件考虑是否要进行数据倒置
        %             % out.depth = -flip(out.depth);
        %             % out.t = flip(out.t, 3);
        %             % out.s = flip(out.s, 3);
        %             % [out.thetao] = dataDepthInterp(out.longitude,out.latitude, out.thetao, out.depth, 0:10:1000);
        %             % [out.so] = dataDepthInterp(out.longitude,out.latitude, out.so, out.depth, 0:10:1000);
        %             % out.depth = 0:10:1000;
        %             %%%%%%%%%%%%%%%%%%
        %             % [obj.pt] = GDC.potiential_temp_measure_from_temp_sail(out.lon, out.lat, out.depth, out.t, out.s);
        %             % [out.pt] = GDC.potiential_temp_measure_from_temp_sail(out.lon, out.lat, out.depth, out.t, out.s);
        %             outArray(i) = out;
        %             % surface = out.thetao(:,:,1);
        %             Sum_t = Sum_t + out.thetao;
        %             Sum_s = Sum_s + out.so;
        % 
        %             t_effictive = ones(size(out.thetao));
        %             t_effictive(isnan(out.thetao)) = 0;
        %             Sum_t_effictive = Sum_t_effictive + t_effictive;
        % 
        % 
        %             s_effictive = ones(size(out.so));
        %             s_effictive(isnan(out.so)) = 0;
        %             Sum_s_effictive = Sum_s_effictive + s_effictive;
        % 
        %         end
        %         mean_pt = Sum_t./Sum_t_effictive;
        %         mean_s = Sum_s./Sum_s_effictive;
        % 
        %         for i = 1:length(filelist)
        % 
        %             date = datetime(str2double(filelist(i).name),tensorDate.Month,tensorDate.Day)
        %             try
        % 
        %                 GDC.GridDataSelect(date);
        %             catch
        %                 warning([datestr(date,'yyyymmdd') '未找到文件']);
        %                 continue;
        %             end
        %             out = outArray(i);
        %             lon = out.longitude;
        %             lat = out.latitude;
        %             depth = out.depth;
        % 
        %             pt = out.thetao;
        %             s = out.so;
        %             ssta = pt - mean_pt;
        %             sssa = s - mean_s;
        % 
        %             filename = [datestr(date,'yyyymmdd') '.mat'];
        %             save(['./' filelist(i).name '/' filename],'lon','lat','depth','pt','s','ssta','sssa');
        %         end
        % 
        % 
        %         if ~exist('./每天气候态数据/', 'dir')
        %             mkdir('./每天气候态数据/');
        %         end
        % 
        %         save(['./每天气候态数据/' num2str(days(tensorDate- startDate) + 1) '.mat'], 'lon','lat','depth','mean_pt','mean_s');
        %         disp('数据保存完成');
        %         tensorDate = tensorDate + days(1);
        %     end
        %     function [vq] = dataDepthInterp(lon, lat, data, old_depth, new_depth)
        %         F = griddedInterpolant({lon, lat, old_depth}, data);
        %         vq = F({lon, lat, new_depth});
        %     end
        % 
        % end

    end
    methods(Static)
        function [lon, lat, depth, pt, s, rho, dh] = parload(file)
            load(file, 'lon','lat','depth','pt','s','rho','dh'); 
        end
        function parsave(filename, lon, lat, depth, pt, s, rho, dh)
            save(filename,'lon','lat','depth','pt','s','rho','dh');
        end

        function parsave2(filename, lon, lat, depth, ssta, sssa, rhoa, dha)
            save(filename,'lon','lat','depth','ssta','sssa','rhoa','dha');
        end

        % 对group进行温盐异常的计算
        function group_measure_and_save(groupfilelist)
            Sum_s = nan(301,301,36);
            Sum_s_effictive = 0;
            Sum_pt = nan(301,301,36);
            Sum_pt_effictive = 0;
            Sum_rho = nan(301,301,36);
            Sum_rho_effictive = 0;
            Sum_dh = nan(301,301,36);
            Sum_dh_effictive = 0;
            for i = 1:length(groupfilelist)
                file = [groupfilelist(i).folder '\' groupfilelist(i).name];
                if ~exist(['./' groupfilelist(i).folder(end-3:end) '/'], 'dir')
                    mkdir(['./' groupfilelist(i).folder(end-3:end) '/']);
                end

                try
                    load(file,'pt','s','rho','dh');

                catch
                    warning([file ' 未找到文件']);
                    ncfilePath = 'F:\ftp\3D sst\Global Ocean Physics Reanalysis\';
                    ncGDC = GridDatasetClass(ncfilePath);
                    [~, ~, ~, pt, s, rho, dh] = ncGDC.rho_dh_measure(i);
                    disp([file ' 已修正']); 
                end 
                [Sum_s, Sum_s_effictive] = sum_measure(s, Sum_s, Sum_s_effictive);
                [Sum_pt, Sum_pt_effictive] = sum_measure(pt, Sum_pt, Sum_pt_effictive);
                [Sum_rho, Sum_rho_effictive] = sum_measure(rho, Sum_rho, Sum_rho_effictive);
                [Sum_dh, Sum_dh_effictive] = sum_measure(dh, Sum_dh, Sum_dh_effictive);
                clear pt s rho dh;
            end

            mean_pt = Sum_pt./Sum_pt_effictive;
            mean_s = Sum_s./Sum_s_effictive;
            mean_rho = Sum_rho./Sum_rho_effictive;
            mean_dh = Sum_dh./Sum_dh_effictive;
            disp('开始保存mat');


            for i = 1:length(groupfilelist)
                file = [groupfilelist(i).folder '\' groupfilelist(i).name];
                try                 
                    [lon, lat, depth, pt, s, rho, dh] = GridDatasetClass.parload(file);       
                catch
                    warning([file ' 未找到文件']);
                    ncfilePath = 'F:\ftp\3D sst\Global Ocean Physics Reanalysis\';
                    ncGDC = GridDatasetClass(ncfilePath);
                    [lon, lat, depth, pt, s, rho, dh] = ncGDC.rho_dh_measure(i);
                    disp([file ' 已修正']); 
                end 
                ssta = pt - mean_pt;
                sssa = s - mean_s;
                rhoa = rho - mean_rho;
                dha = dh - mean_dh;
                % surface = ssta(:,:,23);
                savefilename = ['./' groupfilelist(i).folder(end-3:end) '/' groupfilelist(i).name(1:end-4) '.mat'];

                GridDatasetClass.parsave2(savefilename, lon, lat, depth, ssta, sssa, rhoa, dha);
                disp([savefilename '保存完毕']);
            end
            function [Sum, Sum_effictive] = sum_measure(data, Sum, Sum_effictive)
                % 有效数据计数
                data_effictive = ones(size(data));
                data_effictive(isnan(data)) = 0;
                Sum_effictive = Sum_effictive + data_effictive;
                
               

                % 进行nan处理,把要累加的nan格子变成0
                Sum(~isnan(data)&isnan(Sum)) = 0;
                data(isnan(data)&~isnan(Sum)) = 0;

                % 进行累加
                Sum = Sum + data;
            end


        end
        function [PT,SP,rho,dh] = potiential_temp_measure_from_temp_sail(lon, lat, depth, ptemp, sail)
            [lon,lat,depth] = meshgrid(lon,lat,depth);
            lon = permute(lon,[2 1 3]);
            lat = permute(lat,[2 1 3]);
            depth = permute(depth,[2 1 3]);
            SP = sail;
            PT = ptemp;
            z1 = gsw_z_from_depth(depth);
            p = gsw_p_from_z(z1,lat);
            %% 温度和盐度变换
            [SA, ~] = gsw_SA_from_SP(SP,p,lon,lat);
            %% 位势密度和位势涡度
            CT = gsw_CT_from_pt(SA,PT);
            rho = gsw_rho(SA,CT,p);
            [m,n,z] = size(rho);


            %% 计算动态高度前需要将三维矩阵转化为二维且表现为海洋深度*剖面序列号
            SA_reshape =  reshape(SA,m*n,z)';
            CT_reshape = reshape(CT,m*n,z)';
            p = reshape(p,m*n,z)';

            %% 剔除nan的序列号，参与计算
            idx = ~isnan(SA_reshape(end,:));
            SA_reshape = SA_reshape(:,idx);
            CT_reshape = CT_reshape(:,idx);
            p = p(:,idx);
            if isempty(p)
                dh_effictive = [];
            else
                dh_effictive = gsw_geo_strf_dyn_height(SA_reshape,CT_reshape,p,1000);            
            end
            %% 再恢复剔除序列号，进行矩阵reshape
            dh_T = nan(z,m*n);
            count = 1;
            for i = 1:m*n
                if idx(i)==1
                    dh_T(:,i) = dh_effictive(:,count);
                    count = count + 1;
                end
            end

            dh = reshape(dh_T',m,n,z);

        end

        function [rd_result] = ncFileRead(file, varargin)
            % fileInfo = ncinfo(file);
            % nc_variables = {fileInfo.Variables.Name};
            % for i = 1:length(varargin)
            % 
            %     if ismember(varargin{i},nc_variables)
            %         rd_result.(varargin{i}) = squeeze(ncread(file, varargin{i}));
            %     end
            % end
            % 做判断过于浪费时间
            for i = 1:length(varargin)
                rd_result.(varargin{i}) = squeeze(ncread(file, varargin{i}));
            end
        end


    end
end

